home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / bessel.i < prev    next >
Text File  |  1995-07-26  |  10KB  |  427 lines
  1. /*
  2.    BESSEL.I
  3.    A few Bessel functions.
  4.  
  5.    $Id: bessel.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* Taken from Numerical Recipes.  */
  11.  
  12. /* ------------------------------------------------------------------------ */
  13.  
  14. func bessj0(x)
  15. /* DOCUMENT bessj0(x)
  16.      returns Bessel function J0 at points X.
  17.    SEE ALSO: bessj
  18.  */
  19. {
  20.   ax= abs(x);
  21.   small= (ax<8.0);
  22.   list= where(small);
  23.   if (numberof(list)) {
  24.     xx= x(list);
  25.     y= xx*xx;
  26.     s= poly(y, 57568490574.0, -13362590354.0, 651619640.7, -11214424.18,
  27.         77392.33017, -184.9052456) /
  28.        poly(y, 57568490411.0, 1029532985.0, 9494680.718, 59272.64853,
  29.         267.8532712, 1.0);
  30.   }
  31.   list= where(!small);
  32.   if (numberof(list)) {
  33.     x= x(list);
  34.     ax= abs(x);
  35.     z= 8.0/ax;
  36.     y= z*z;
  37.     xx= ax-0.785398164;  /* pi/4, rounded incorrectly */
  38.     l= sqrt(0.636619772/ax) *
  39.        (cos(xx)*poly(y, 1.0, -0.1098628627e-2,
  40.              0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6) -
  41.     sin(xx)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3,
  42.                -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7));
  43.   }
  44.   return merge(s, l, small);
  45. }
  46.  
  47. func bessj1(x)
  48. /* DOCUMENT bessj1(x)
  49.      returns Bessel function J1 at points X.
  50.    SEE ALSO: bessj
  51.  */
  52. {
  53.   ax= abs(x);
  54.   small= (ax<8.0);
  55.   list= where(small);
  56.   if (numberof(list)) {
  57.     xx= x(list);
  58.     y= xx*xx;
  59.     s= xx * poly(y, 72362614232.0, -7895059235.0, 242396853.1, -2972611.439,
  60.          15704.48260, -30.16036606) /
  61.        poly(y, 144725228442.0, 2300535178.0, 18583304.74, 99447.43394,
  62.         376.9991397, 1.0);
  63.   }
  64.   list= where(!small);
  65.   if (numberof(list)) {
  66.     x= x(list);
  67.     ax= abs(x);
  68.     z= 8.0/ax;
  69.     y= z*z;
  70.     xx= ax-2.356194491;  /* 3*pi/4 */
  71.     l= sign(x) * sqrt(0.636619772/ax) *
  72.        (cos(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4,
  73.              0.2457520174e-5, -0.240337019e-6) -
  74.     sin(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5,
  75.                -0.88228987e-6, 0.105787412e-6));
  76.   }
  77.   return merge(s, l, small);
  78. }
  79.  
  80. func bessj(n, x)
  81. /* DOCUMENT bessj(n, x)
  82.      returns Bessel function Jn of order N at points X.  N must be scalar.
  83.    SEE ALSO: bessy, bessi, bessk, bessj0, bessj1
  84.  */
  85. {
  86.   if (n>1) {
  87.     ax= abs(x);
  88.     big= (ax > n);
  89.     list= where(big);
  90.     if (numberof(list)) {
  91.       /* upward recurrence */
  92.       ax= abs(x(list));
  93.       tox= 2.0/ax;
  94.       bjm= bessj0(ax);
  95.       bj= bessj1(ax);
  96.       for (i=1 ; i<n ; i++) {
  97.     bjp= i*tox*bj-bjm;
  98.     bjm= bj;
  99.     bj= bjp;
  100.       }
  101.     }
  102.     list= where(!big);
  103.     if (numberof(list)) {
  104.       ax= abs(x(list));
  105.       zero= (ax==0.0);
  106.       list= where(zero);
  107.       if (numberof(list)) {
  108.     bj0= ax(list);  /* == 0.0 */
  109.       }
  110.       list= where(!zero);
  111.       if (numberof(list)) {
  112.     /* downward recurrence */
  113.     ax= ax(list);
  114.     tox= 2.0/ax;
  115.     m= 2*((n+long(sqrt(bess_acc*n)))/2);
  116.     jsum= 0;
  117.     bjp= ans= add= array(0.0, numberof(ax));
  118.     bj1= array(1.0, numberof(ax));
  119.     for (i=m ; i>0 ; i--) {
  120.       bjm= i*tox*bj1-bjp;
  121.       bjp= bj1;
  122.       bj1= bjm;
  123.       renorm= (abs(bj1)>bess_big);
  124.       list= where(renorm);
  125.       if (numberof(list)) {
  126.         bj1(list)/= bess_big;
  127.         bjp(list)/= bess_big;
  128.         ans(list)/= bess_big;
  129.         add(list)/= bess_big;
  130.       }
  131.       if (jsum) add+= bj1;
  132.       jsum= !jsum;
  133.       if (i==n) ans= bjp;
  134.     }
  135.     bj1= ans/(2.0*add-bj1);
  136.       }
  137.       bj1= merge(bj0, bj1, zero);
  138.     }
  139.     bj= merge(bj, bj1, big);
  140.     if (n%2) bj*= sign(x);
  141.     return bj;
  142.   } else if (n==1) {
  143.     return bessj1(x);
  144.   } else if (!n) {
  145.     return bessj0(x);
  146.   }
  147. }
  148.  
  149. /* ------------------------------------------------------------------------ */
  150.  
  151. func bessy0(x)
  152. /* DOCUMENT bessy0(x)
  153.      returns Bessel function Y0 at points X.
  154.    SEE ALSO: bessy
  155.  */
  156. {
  157.   ax= abs(x);
  158.   small= (ax<8.0);
  159.   list= where(small);
  160.   if (numberof(list)) {
  161.     xx= x(list);
  162.     y= xx*xx;
  163.     s= poly(y, -2957821389.0, 7062834065.0, -512359803.6, 10879881.29,
  164.         -86327.92757, 228.4622733) /
  165.        poly(y, 40076544269.0, 745249964.8, 7189466.438, 47447.26470,
  166.         226.1030244, 1.0) + 0.636619772*bessj0(x)*log(x);
  167.   }
  168.   list= where(!small);
  169.   if (numberof(list)) {
  170.     x= x(list);
  171.     ax= abs(x);
  172.     z= 8.0/ax;
  173.     y= z*z;
  174.     xx= ax-0.785398164;  /* pi/4, rounded incorrectly */
  175.     l= sqrt(0.636619772/ax) *
  176.        (sin(xx)*poly(y, 1.0, -0.1098628627e-2, 0.2734510407e-4,
  177.              -0.2073370639e-5, 0.2093887211e-6) -
  178.     cos(xx)*z*poly(y, -0.1562499995e-1, 0.1430488765e-3,
  179.                -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7));
  180.   }
  181.   return merge(s, l, small);
  182. }
  183.  
  184. func bessy1(x)
  185. /* DOCUMENT bessy1(x)
  186.      returns Bessel function Y1 at points X.
  187.    SEE ALSO: bessy
  188.  */
  189. {
  190.   ax= abs(x);
  191.   small= (ax<8.0);
  192.   list= where(small);
  193.   if (numberof(list)) {
  194.     xx= x(list);
  195.     y= xx*xx;
  196.     s= xx * poly(y, -0.4900604943e13, 0.1275274390e13, -0.5153438139e11,
  197.          0.7349264551e9, -0.4237922726e7, 0.8511937935e4) /
  198.        poly(y, 0.2499580570e14, 0.4244419664e12, 0.3733650367e10,
  199.         0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1.0) +
  200.        0.636619772*(bessj1(x)*log(x)-1.0/x);
  201.   }
  202.   list= where(!small);
  203.   if (numberof(list)) {
  204.     x= x(list);
  205.     ax= abs(x);
  206.     z= 8.0/ax;
  207.     y= z*z;
  208.     xx= ax-2.356194491;  /* 3*pi/4 */
  209.     l= sqrt(0.636619772/x) *
  210.        (sin(xx)*poly(y, 1.0, 0.183105e-2, -0.3516396496e-4,
  211.              0.2457520174e-5, -0.240337019e-6) +
  212.     cos(xx)*z*poly(y, 0.04687499995, -0.2002690873e-3, 0.8449199096e-5,
  213.                -0.88228987e-6, 0.105787412e-6));
  214.   }
  215.   return merge(s, l, small);
  216. }
  217.  
  218. func bessy(n, x)
  219. /* DOCUMENT bessy(n, x)
  220.      returns Bessel function Yn of order N at points X.  N must be scalar.
  221.    SEE ALSO: bessj, bessi, bessk, bessy0, bessy1
  222.  */
  223. {
  224.   if (n>1) {
  225.     /* upward recurrence */
  226.     tox= 2.0/x;
  227.     bym= bessy0(x);
  228.     by= bessy1(x);
  229.     for (i=1 ; i<n ; i++) {
  230.       byp= i*tox*by-bym;
  231.       bym= by;
  232.       by= byp;
  233.     }
  234.     return by;
  235.   } else if (n==1) {
  236.     return bessy1(x);
  237.   } else if (!n) {
  238.     return bessy0(x);
  239.   }
  240. }
  241.  
  242. /* ------------------------------------------------------------------------ */
  243.  
  244. func bessi0(x)
  245. /* DOCUMENT bessi0(x)
  246.      returns Bessel function I0 at points X.
  247.    SEE ALSO: bessi
  248.  */
  249. {
  250.   ax= abs(x);
  251.   small= (ax<3.75);
  252.   list= where(small);
  253.   if (numberof(list)) {
  254.     xx= x(list)/3.75;
  255.     y= xx*xx;
  256.     s= poly(y, 1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732,
  257.         0.360768e-1, 0.45813e-2);
  258.   }
  259.   list= where(!small);
  260.   if (numberof(list)) {
  261.     x= x(list);
  262.     ax= abs(x);
  263.     y= 3.75/ax;
  264.     l= (exp(ax)/sqrt(ax)) * poly(y, 0.39894228, 0.1328592e-1, 0.225319e-2,
  265.                  -0.157565e-2, 0.916281e-2, -0.2057706e-1,
  266.                  0.2635537e-1, -0.1647633e-1, 0.392377e-2);
  267.   }
  268.   return merge(s, l, small);
  269. }
  270.  
  271. func bessi1(x)
  272. /* DOCUMENT bessi1(x)
  273.      returns Bessel function I1 at points X.
  274.    SEE ALSO: bessi
  275.  */
  276. {
  277.   ax= abs(x);
  278.   small= (ax<3.75);
  279.   list= where(small);
  280.   if (numberof(list)) {
  281.     xx= x(list);
  282.     y= xx/3.75;
  283.     y*= y;
  284.     s= abs(xx) * poly(y, 0.5, 0.87890594, 0.51498869, 0.15084934,
  285.               0.2658733e-1, 0.301532e-2, 0.32411e-3);
  286.   }
  287.   list= where(!small);
  288.   if (numberof(list)) {
  289.     ax= ax(list);
  290.     y= 3.75/ax;
  291.     l= (exp(ax)/sqrt(ax)) *
  292.        poly(y, 0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2,
  293.         -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1,
  294.         -0.420059e-2);
  295.   }
  296.   return sign(x) * merge(s, l, small);
  297. }
  298.  
  299. func bessi(n, x)
  300. /* DOCUMENT bessi(n, x)
  301.      returns Bessel function In of order N at points X.  N must be scalar.
  302.    SEE ALSO: bessk, bessj, bessy, bessi0, bessi1
  303.  */
  304. {
  305.   if (n>1) {
  306.     zero= (x==0.0);
  307.     list= where(zero);
  308.     if (numberof(list)) {
  309.       bi0= x(list);  /* == 0.0 */
  310.     }
  311.     list= where(!zero);
  312.     if (numberof(list)) {
  313.       /* downward recurrence */
  314.       x= x(list);
  315.       ax= abs(x);
  316.       tox= 2.0/ax;
  317.       m= 2*(n+long(sqrt(bess_acc*n)));
  318.       bip= ans= array(0.0, numberof(ax));
  319.       bi= array(1.0, numberof(ax));
  320.       for (i=m ; i>0 ; i--) {
  321.     bim= i*tox*bi+bip;
  322.     bip= bi;
  323.     bi= bim;
  324.     list= where(abs(bi) > bess_big);
  325.     if (numberof(list)) {
  326.       ans(list)/= bess_big;
  327.       bi(list)/= bess_big;
  328.       bip(list)/= bess_big;
  329.     }
  330.     if (i==n) ans= bip;
  331.       }
  332.       bi= ans*bessi0(x)/bi;
  333.       if (n%2) bi*= sign(x);
  334.     }
  335.     return merge(bi0, bi, zero);
  336.   } else if (n==1) {
  337.     return bessi1(x);
  338.   } else if (!n) {
  339.     return bessi0(x);
  340.   }
  341. }
  342.  
  343. /* ------------------------------------------------------------------------ */
  344.  
  345. func bessk0(x)
  346. /* DOCUMENT bessk0(x)
  347.      returns Bessel function K0 at points X.
  348.    SEE ALSO: bessk
  349.  */
  350. {
  351.   small= (x<=2.0);
  352.   list= where(small);
  353.   if (numberof(list)) {
  354.     xx= x(list);
  355.     y= xx*xx/4.0;
  356.     s= (-log(xx/2.0)*bessi0(xx)) +
  357.        poly(y, -0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2,
  358.         0.10750e-3, 0.74e-5);
  359.   }
  360.   list= where(!small);
  361.   if (numberof(list)) {
  362.     x= x(list);
  363.     y= 2.0/x;
  364.     l= (exp(-x)/sqrt(x)) *
  365.        poly(y, 1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1,
  366.         0.587872e-2, -0.251540e-2, 0.53208e-3);
  367.   }
  368.   return merge(s, l, small);
  369. }
  370.  
  371. func bessk1(x)
  372. /* DOCUMENT bessk1(x)
  373.      returns Bessel function K1 at points X.
  374.    SEE ALSO: bessk
  375.  */
  376. {
  377.   small= (x<=2.0);
  378.   list= where(small);
  379.   if (numberof(list)) {
  380.     xx= x(list);
  381.     y= xx*xx/4.0;
  382.     s= (log(xx/2.0)*bessi1(xx)) +
  383.        (1.0/xx) * poly(y, 1.0, 0.15443144, -0.67278579, -0.18156897,
  384.                -0.1919402e-1, -0.110404e-2, -0.4686e-4);
  385.   }
  386.   list= where(!small);
  387.   if (numberof(list)) {
  388.     x= x(list);
  389.     y= 2.0/x;
  390.     l= (exp(-x)/sqrt(x)) *
  391.        poly(y, 1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1,
  392.         -0.780353e-2, 0.325614e-2, -0.68245e-3);
  393.   }
  394.   return merge(s, l, small);
  395. }
  396.  
  397. func bessk(n, x)
  398. /* DOCUMENT bessk(n, x)
  399.      returns Bessel function Kn of order N at points X.  N must be scalar.
  400.    SEE ALSO: bessi, bessj, bessy, bessi0, bessi1
  401.  */
  402. {
  403.   if (n>1) {
  404.     /* upward recurrence */
  405.     tox= 2.0/x;
  406.     bkm= bessk0(x);
  407.     bk= bessk1(x);
  408.     for (i=1 ; i<n ; i++) {
  409.       bkp= i*tox*bk+bkm;
  410.       bkm= bk;
  411.       bk= bkp;
  412.     }
  413.     return bk;
  414.   } else if (n==1) {
  415.     return bessk1(x);
  416.   } else if (!n) {
  417.     return bessk0(x);
  418.   }
  419. }
  420.  
  421. /* ------------------------------------------------------------------------ */
  422.  
  423. bess_acc= 40.0;
  424. bess_big= 1.e10;
  425.  
  426. /* ------------------------------------------------------------------------ */
  427.